#noteB1_FindMultiModalNeighbors()



(1) FindMultiModalNeighbors() 构建带权重的最近邻图(WNN)

$ find . | grep -v "dustbin/" | grep -v "notes/" | grep "\.R$" | xargs grep -in "FindMultiModalNeighbors" --color=auto 2>/dev/null
./seurat-4.4.0/R/clustering.R:45:FindMultiModalNeighbors  <- function(

#' Construct weighted nearest neighbor graph
#'
#' This function will construct a weighted nearest neighbor (WNN) graph. 
#' 本函数构建权重最近邻(WNN)图。
#'
#' For each cell, we identify the nearest neighbors based on a weighted combination of two modalities. 
#' 对每个细胞，基于两个权重的权重组合，识别出最近邻。
#'
#' Takes as input two dimensional reductions, one computed for each modality.
#' 输入是两个降维结果，对应两个模态。
#' 
#' Other parameters are listed for debugging, but can be left as default values.
#' 其他参数列举如下，用于排错，但可以不管直接使用默认值。
#'
#' @param object A Seurat object
#' @param reduction.list A list of two dimensional reductions, one for each of
#' the modalities to be integrated | 两个降维列表，分别对应用于整合的两个模态
#' @param dims.list A list containing the dimensions for each reduction to use | 每个降维使用的维度
#' @param k.nn the number of multimodal neighbors to compute. 20 by default | 找多模态最近邻的个数，默认20
#' @param l2.norm Perform L2 normalization on the cell embeddings after
#' dimensional reduction. TRUE by default. | 对降维后的细胞嵌入做L2标准化，默认T
#' @param knn.graph.name Multimodal knn graph name | 多模态knn图名字
#' @param snn.graph.name Multimodal snn graph name | 多模态snn图名字
#' @param weighted.nn.name Multimodal neighbor object name | 多模态邻居对象名字
#' @param modality.weight.name Variable name to store modality weight in object
#' meta data | 在meta.data中保存模态权重的变量名
#' @param knn.range The number of approximate neighbors to compute | 要计算的近似邻居的数量
#' @param prune.SNN Cutoff not to discard edge in SNN graph | SNN图中丢掉边的截断值
#' @param sd.scale  The scaling factor for kernel width. 1 by default | 和函数宽度的缩放因子，默认1
#' @param cross.contant.list Constant used to avoid divide-by-zero errors. 1e-4
#' by default | 常量，用于防止除以0错误。 默认 1e-4
#' @param smooth Smoothing modality score across each individual modality
#' neighbors. FALSE by default | 平滑每个模态内相邻邻居的模态得分。默认F
#' @param return.intermediate Store intermediate results in misc | 将中间结果存储在misc中
#' @param modality.weight A \code{\link{ModalityWeights}} object generated by
#' \code{FindModalityWeights}  | FindModalityWeights 产生的 ModalityWeights 对象
#' @param verbose Print progress bars and output | 打印进度条 和 输出
#'
#' @return Seurat object containing a nearest-neighbor object, KNN graph, and
#' SNN graph - each based on a weighted combination of modalities. | 返回Seurat对象，包含一个最近邻对象，knn图，snn图，都是模态的权重组合。
#' @concept clustering
#' @export
#'
FindMultiModalNeighbors  <- function(
  object,
  reduction.list,
  dims.list,
  k.nn = 20,
  l2.norm = TRUE,
  knn.graph.name = "wknn",
  snn.graph.name = "wsnn",
  weighted.nn.name = "weighted.nn",
  modality.weight.name = NULL,
  knn.range = 200,
  prune.SNN = 1/15,
  sd.scale = 1,
  cross.contant.list = NULL,
  smooth = FALSE,
  return.intermediate = FALSE,
  modality.weight = NULL,
  verbose = TRUE
) {
  #(A1) 如果为空，设置默认值为 1e-4
  cross.contant.list <- cross.contant.list %||%
    as.list(x = rep(x = 1e-4, times = length(x = reduction.list)))
    # reduction.list 难道还能不是2? 那就是多个模态组合了？

  #(A2) 如果modality.weight为空(默认空)
  if (is.null(x = modality.weight)) {
    if (verbose) {
      # 消息：计算细胞特异的模态权重
      message("Calculating cell-specific modality weights")
    }
    # (B1) 调用函数：计算模态权重，详见下文(2)
    modality.weight <- FindModalityWeights(
      object = object,
      reduction.list = reduction.list,
      dims.list = dims.list,
      k.nn = k.nn,
      sd.scale = sd.scale,
      l2.norm = l2.norm,
      cross.contant.list = cross.contant.list,
      smooth = smooth,
      verbose = verbose
    )
  }


  #(A3) 权重的名字：使用的assay加上.weight后缀
  # 实例输入 modality.weight.name = "RNA.weight"，貌似不够用
  modality.weight.name <- modality.weight.name %||% paste0(modality.weight@modality.assay, ".weight")

  #(A4) 获取 FindModalityWeights 返回中使用的 assay 名字
  modality.assay <- slot(object = modality.weight, name = "modality.assay")
  #  pca  apca 
  #"RNA" "ADT"


  #(A5) 如果名词长度不一致，则警告，并使用默认的设置： 权重返回的 assay 名字，加上 .weight 后缀
  if (length(modality.weight.name) != length(reduction.list)) {
    warning("The number of provided modality.weight.name is not equal to the number of modalities. ",
            paste(paste0(modality.assay, ".weight"), collapse = " "),
            " are used to store the modality weights" )
    modality.weight.name <- paste0(modality.assay, ".weight")
  }


  #(A5) 第一个实验，这个变量在(A17) 采用
  first.assay <- modality.assay[1]

  #(A6) 多模态最近邻 20个：返回cell 下标 和 多模态最近邻距离（欧氏距离的指数函数 * 模态权重）
  weighted.nn <- MultiModalNN( #函数见下文(12)
    object = object,
    k.nn = k.nn,  #20个最近邻
    modality.weight = modality.weight, #模态权重
    knn.range = knn.range, # 每个模态的最近邻数量初始值200个细胞，合并后200-400个，加权距离后排序，取前20个
    verbose = verbose #默认T
  )
  #(A7) 下标
  select_nn <- Indices(object = weighted.nn)
  #(A8) 多模态距离
  select_nn_dist <- Distances(object = weighted.nn)





  #(A9) 计算 多模态KNN图
  # compute KNN graph
  if (verbose) {
    message("Constructing multimodal KNN graph")
  }
  #(A10) 细胞下标，转置，变为一行
  j <- as.numeric(x = t(x = select_nn ))
  # 转置后：列为细胞，及该细胞最近的细胞下标
  # 变为一行：
  # > head(j) #[1] 2014 4815 5515 4135 6050 6473
  # > length(j) #[1] 160720

  #(A11) 细胞的下标编号：序号-1，整除 20，然后加1
  i <- ((1:length(x = j)) - 1) %/% k.nn + 1
  # > head(i)
  #[1] 1 1 1 1 1 1
  #> unique(i) |> length()
  #[1] 8036

  #(A12) 构建大型稀疏矩阵：i行为细胞，j列表示每个细胞下标，值为1
  nn.matrix <- sparseMatrix(
    i = i,
    j = j,
    x = 1,
    dims = c(ncol(x = object), ncol(x = object))
  )

  #(A13)对角线填充1
  diag(x = nn.matrix) <- 1

  #(A14) 行列名都设置为 cell id
  rownames(x = nn.matrix) <-  colnames(x = nn.matrix) <- colnames(x = object)

  #(A15) 这个公式的目的是什么？
  # 我认为这一步不是必须的，跳过会怎样？本步骤可能导致有些cell的最近邻超过20个。 //todo
  # 检查之前的 本系列 源码解析21-2(3) FindNeighbors(A9) 没有这个变换。 直接在稀疏矩阵变图了。
  # 更关键的是：接着计算SNN用的是select_nn，而不是这个经过此公式变换的图！
  # 那这个图记录到Seurat对象中有啥用呢？bm2@graphs$wknn 貌似没用。为何不记录原始计算SNN的KNN图？
  #
  # 分析该公式做了啥：
  # 假设A是元素为0或者1的矩阵，且对角线为1，则经过 A + A^T - A^T.A，得到什么？
  # 结果就是：得到对称矩阵，且对称位置aij, aji只要有一个为1，结果为1。同为0，结果才是0。
  # 证明过程：每一项都是 Cji =Aji + Aij − (Aij ⋅ Aji)，
  # 矩阵的值只有1或者0，那么只有4种类情况：都是0，结果是0。其他3种情况都是1。
  nn.matrix <- nn.matrix + t(x = nn.matrix) - t(x = nn.matrix) * nn.matrix

  #(A16) 转为图类型
  nn.matrix <- as.Graph(x = nn.matrix)
  #(A17)定义slot assay.used，为 modality.assay[1]
  slot(object = nn.matrix, name = "assay.used") <- first.assay

  #(A18) 写入图到obj中: 图类型会保存到 @graphs 中
  object[[knn.graph.name]] <- nn.matrix
  # > bm2@graphs |> names()
  # [1] "wknn" "wsnn"





  #(A19) 消息计算 SNN 图
  # compute SNN graph
  if (verbose) {
    message("Constructing multimodal SNN graph")
  }

  #(A20) 输入每个细胞的20个 多模态最近邻 下标，计算Jacard相似性，小于prune.SNN=1/15的记为0
  snn.matrix <- ComputeSNN(nn_ranked = select_nn, prune = prune.SNN) 
  #ComputeSNN()见本系列 源码解析21-2(13)

  #(A21) 行列名都是cell id
  rownames(x = snn.matrix) <- colnames(x = snn.matrix) <- Cells(x = object)

  #(A22) 转为图
  snn.matrix <- as.Graph(x = snn.matrix )
  # 定义slot assay.used，为 modality.assay[1]
  slot(object = snn.matrix, name = "assay.used") <- first.assay

  #(A23) 写入图到obj中: 图类型会保存到 @graphs 中
  object[[snn.graph.name]] <- snn.matrix



  #(A24) 写入图到obj中: Neighbor 类型会保存到 @neighbors 中
  # add neighbors and modality weights
  object[[weighted.nn.name]] <- weighted.nn
  # > bm2@neighbors |> names()
  # [1] "weighted.nn"


  #(A25) 把模态权重保存到 meta.data 中的两列：RNA.weight, ADT.weight
  for (m in 1:length(x = modality.weight.name)) {
    object[[modality.weight.name[[m]]]] <- slot(
      object = modality.weight,
      name = "modality.weight.list"
    )[[m]]
  }
  #> slot(
  #     object = modality.weight_2,
  #     name = "modality.weight.list"
  # )[[1]] |> head(n=3)
  # AAACCTGAGACTCGGA-1 AAACCTGAGAGTACAT-1 AAACCTGAGATATACG-1 
  #          0.5325040          0.4924699          0.5799385 


  #(A26) 记录日志
  # add command log
  #(B1) 获取 模态权重计算中的 command slot
  modality.weight.command <- slot(object = modality.weight, name = "command")
  #(B2) 设置 assay.used 为第一个模态
  slot(object = modality.weight.command, name = "assay.used") <- first.assay
  #(B3) 获取 name： @ name       : chr "Seurat:::FindModalityWeights"
  modality.weight.command.name <- slot(object = modality.weight.command, name = "name")
  #(B4) 保存到对象中，该对象要返回的
  object[[modality.weight.command.name]] <- modality.weight.command
  
  #(B5) ？？ 不确定，函数外还不能调用。猜测：获取同名记录，更新后插入回去
  # Error: 'LogSeuratCommand' cannot be called at the top level
  command <- LogSeuratCommand(object = object, return.command = TRUE)
  slot(object = command, name = "params")$modality.weight <- NULL
  slot(object = command, name = "assay.used") <- first.assay
  command.name <- slot(object = command, name = "name")
  object[[command.name]] <- command

  #(A27) 如果返回中间结果，放在Misc中。默认F
  if (return.intermediate) {
    Misc(object = object, slot = "modality.weight") <- modality.weight
  }
  return (object)
}









(2) FindModalityWeights() 计算模态的权重
$ find . | grep -v "dustbin/" | grep -v "notes/" | xargs grep -in "FindModalityWeights" --color=auto 2>/dev/null
./seurat-4.4.0/R/clustering.R:1028:FindModalityWeights  <- function(

* 这个函数很复杂，涉及好几个函数，最基层又是c++函数。
* 这次对c++函数做了初步测试，见 https://blog.csdn.net/wangjunliang/article/details/143096799
* 尝试阅读3遍
  - 第一遍：逐行阅读注释，尝试仅用部分数据做测试
  - 第二遍：整个函数执行一次，再次按代码块看，从结果倒推、修正第一遍的错误
  - 第三遍：结合其他人的心得体会、论文描述、会议见闻，反思、抽象，给出整体处理过程的描述。


# Calculate modality weights | 计算模块权重，逐行阅读
#
# This function calculates cell-specific modality weights which are used to
# in WNN analysis. | 该函数计算细胞特异的模态权重，用于WNN分析
#' @inheritParams FindMultiModalNeighbors |参数来源于 FindMultiModalNeighbors()
# @param object A Seurat object
# @param snn.far.nn Use SNN farthest neighbors to calculate the kernel width | 使用SNN最远邻居计算核函数宽度
# @param s.nn How many SNN neighbors to use in kernel width | 现有核宽度内，使用的SNN邻居个数
# @param sigma.idx Neighbor index used to calculate kernel width if snn.far.nn = FALSE | 如果snn.far.nn=F使用邻居索引下标计算核宽度
# @importFrom pbapply pblapply | 从pbapply包引用pblapply函数
# @return Returns a \code{ModalityWeights} object that can be used as input to
# \code{\link{FindMultiModalNeighbors}} | 返回值是一个 ModalityWeights 对象，可以作为 FindMultiModalNeighbors 的输入
#
#' @importFrom pbapply pblapply
#
FindModalityWeights  <- function(
  object,
  reduction.list,
  dims.list,
  k.nn = 20,
  snn.far.nn = TRUE,
  s.nn = k.nn,
  prune.SNN = 0,
  l2.norm = TRUE,
  sd.scale = 1,
  query = NULL,
  cross.contant.list = NULL,
  sigma.idx = k.nn,
  smooth = FALSE,
  verbose = TRUE
) {
  # 说明：这几个函数不好理解，都带入实例数据具体化了（RNA+ADT多组学整合）

  #(A1) verbose 默认T，则 my.lapply 使用 pblapply，比普通lapply 多了pb(progress bar进度条)
  my.lapply <- ifelse(
    test = verbose,
    yes = pblapply,
    no = lapply
  )

  #(A2) 默认值 1e-4，说是防止除以0错误: 长度>=2
  cross.contant.list <- cross.contant.list %||%
    as.list(x = rep(x = 1e-4, times = length(x = reduction.list)))
  # > as.list(x = rep(x = 1e-4, times = 2))
#[[1]]
#[1] 1e-04
#
#[[2]]
#[1] 1e-04

  #(A3) 拆开list，变成向量： c("pca", "apca")
  reduction.set <- unlist(x = reduction.list)
  
  #(A4) 连写法！给几个list的键命名
  names(x = reduction.list) <- names(x = dims.list) <-
    names(x = cross.contant.list) <- reduction.set
  
  #(A5) 获取降维的细胞嵌入：每个模态的前若干个PC
  embeddings.list <- lapply(
    X = reduction.list,
    FUN = function(r) Embeddings(object = object, reduction = r)[, dims.list[[r]]]
  )
#> str(embeddings.list)
#List of 2
# $ pca : num [1:8036, 1:30] 0.458 0.316 -0.808 6.787 3.206 ...
#  ..- attr(*, "dimnames")=List of 2
#  .. ..$ : chr [1:8036] "AAACCTGAGACTCGGA-1" "AAACCTGAGAGTACAT-1" "AAACCTGAGATATACG-1" "AAACCTGAGCGAAGGG-1" ...
#  .. ..$ : chr [1:30] "PC_1" "PC_2" "PC_3" "PC_4" ...
# $ apca: num [1:8036, 1:20] -2.606 -3.207 -3.187 3.32 -0.354 ...
#  ..- attr(*, "dimnames")=List of 2
#  .. ..$ : chr [1:8036] "AAACCTGAGACTCGGA-1" "AAACCTGAGAGTACAT-1" "AAACCTGAGATATACG-1" "AAACCTGAGCGAAGGG-1" ...
#  .. ..$ : chr [1:20] "apca_1" "apca_2" "apca_3" "apca_4" ...
#

  #(A6) 默认T，做L2标准化
  if (l2.norm) {
  	#(B1) 对每个矩阵，分别做L2标准化，该函数 L2Norm 见下文(3)
    embeddings.list.norm <- lapply(
      X = embeddings.list,
      # L2Norm() 默认按行，除以每行的平方和的开方，相当于现在每行的cell都在一个单位1球面上
      FUN = function(embeddings) L2Norm(mat = embeddings)
    )
  } else {
  	#(B2) 默认不走这里。
    embeddings.list.norm <- embeddings.list
  }
#实例数据：
#embeddings.list.norm <- lapply(
#  X = embeddings.list,
#  FUN = function(embeddings){
#    print(dim(embeddings))
#    print(embeddings[1:2, 1:3])
#    res=Seurat:::L2Norm(mat = embeddings)
#    print(res[1:2, 1:3])
#    res
#  }
#)
#[1] 8036   30
#                        PC_1       PC_2      PC_3
#AAACCTGAGACTCGGA-1 0.4580905  -8.923786 2.0374376
#AAACCTGAGAGTACAT-1 0.3162329 -10.242548 0.7259737
# 输出是相似的结构
#                         PC_1       PC_2       PC_3
#AAACCTGAGACTCGGA-1 0.03717500 -0.7241838 0.16534230
#AAACCTGAGAGTACAT-1 0.02301086 -0.7453045 0.05282587
#
#[1] 8036   20
#                      apca_1    apca_2    apca_3
#AAACCTGAGACTCGGA-1 -2.605531 -1.501454 -2.837464
#AAACCTGAGAGTACAT-1 -3.207008 -1.267772  4.511293
# 输出是相似的结构
#                       apca_1     apca_2     apca_3
#AAACCTGAGACTCGGA-1 -0.4475578 -0.2579081 -0.4873974
#AAACCTGAGAGTACAT-1 -0.4730865 -0.1870172  0.6654899
# > all( abs(rowSums(embeddings.list.norm$pca**2) - 1) < 1e-10 ) #T


  #(A7) query 默认为null，则定义其为 obj
  if (is.null(x = query)) {
  	# (B1)默认走这: 
    query.embeddings.list.norm <- embeddings.list.norm
    query <- object
  } else {
  	# (B2) 默认不走这
    if (snn.far.nn) {
      stop("query does not support using snn to find distant neighbors")
    }
    query.embeddings.list <- lapply(
      X = reduction.list,
      FUN = function(r) {
        Embeddings(object = query, reduction = r)[, dims.list[[r]]]
      }
    )
    if (l2.norm) {
      query.embeddings.list <- lapply(
        X = query.embeddings.list,
        FUN = function(embeddings) L2Norm(mat = embeddings)
      )
    }
    query.embeddings.list.norm <- query.embeddings.list
  }


  #(A8) 消息：每个模态，找20个最近邻
  if (verbose) {
    message("Finding ", k.nn, " nearest neighbors for each modality.")
  }

  #(A9) 模态内，用指定距离，求每行(obs)的20个最近邻居(Neighbor对象)，给出下标和距离
  nn.list <- my.lapply(
    X = reduction.list,
    FUN = function(r) {
      nn.r <- NNHelper( # NNHelper 该函数的解析见本系列解析21-2(4)
        data = embeddings.list.norm[[r]],
        query = query.embeddings.list.norm[[r]],
        k = max(k.nn, sigma.idx, s.nn),
        method = "annoy",
        metric = "euclidean"
      )
      return(nn.r)
    }
  )

  #(A10) 最近邻对象 赋值 给新变量名
  sigma.nn.list <- nn.list
  
  #(A11) 如果k.nn小，则截断到k.nn个最近邻
  #默认这三个值相等，都等于k.nn：则跳过
  if (sigma.idx > k.nn || s.nn > k.nn) {
    nn.list <- lapply(
      X = nn.list,
      FUN = function(nn){
        slot(object = nn, name = "nn.idx") <- Indices(object = nn)[, 1:k.nn]
        slot(object = nn, name = "nn.dists") <- Distances(object = nn)[, 1:k.nn]
        return(nn)
      }
    )
  }

  #(A12) 获取和最近临的一个细胞的距离
  # Distances() 函数：输入一个 Neighbor 对象，返回啥？返回带cell id行名的 距离矩阵。 该函数详见下文(4)
  # 去第二列，就是除自己之外的最近的邻居的距离
  nearest_dist <- lapply(X = reduction.list, FUN = function(r) Distances(object = nn.list[[r]])[, 2])
  

  #(A13) 预测模态内和模态间的细胞嵌入 【重点】
  # 根据模态内和间的最近邻，计算在原空间的 cell embeding(normPC_表示) 均值
  within_impute <- list()
  cross_impute <- list()
  # Calculating within and cross modality distance
  for (r in reduction.set) { #"pca"  "apca"
  	#(B1) 遍历 模态类型：c("pca", "apca") 创建新名字，模态名 + .norm结尾
    reduction.norm <- paste0(r, ".norm")

    #(B2) 在对象中添加新降维对象: 见本系列 源码解析19-2(2)
    # 产生新的降维类 "pca.norm": 降维必须知道对应的是哪个assay
    object[[ reduction.norm ]] <- CreateDimReducObject(
      embeddings = embeddings.list.norm[[r]], #放到了新对象的 @cell.embeddings 中
      key = paste0("norm", Key(object = object[[r]])), #normPC_
      assay = DefaultAssay(object = object[[r]]) #RNA
    )

    #(B3) 模态内预测： 最近邻的19个细胞（去掉自己）的下标，在自己L2标准化PC空间中求均值，是其在 normPC_x 空间的预测值
    # PredictAssay 函数详见(5)
    within_impute[[r]] <- PredictAssay(
      object = object,
      nn.idx =  Indices(object = nn.list[[r]]), #返回最近邻索引 cell*k 矩阵，行名是cel id。见(6)。是 'pca' 空间内
      reduction = reduction.norm, # "pca.norm" 取20个细胞的嵌入均值时，使用的数据名。'pca' 空间内
      dims = 1:ncol(x = embeddings.list.norm[[r]]), #1:PC个数
      verbose = FALSE,
      return.assay = FALSE
    )
    # > within_impute[[r]][1:2, 1:3]
    #         AAACCTGAGACTCGGA-1 AAACCTGAGAGTACAT-1 AAACCTGAGATATACG-1
    #normPC_1         0.01897473         -0.0372505         -0.0606472
    #normPC_2        -0.82485107         -0.7496652         -0.7927484


    #(B4) 跨模态预测：除了本模态外的模态。
    # 问题： 是否支持 >=3 个模态的整合？ //todo
    cross.modality <- setdiff(x = reduction.set, y = r) #"cpca"

    #(B5) 模态之间预测： 根据其他模态 获取 最近邻的19个细胞（去掉自己），在自己原L2标准化PC空间求均值，是其在原 normPC_x 空间的预测值
    cross_impute[[r]] <- lapply(X = cross.modality, FUN = function(r2) {
      PredictAssay(
        object = object,
        nn.idx = Indices(object = nn.list[[r2]]),  #返回最近邻索引 cell*k 矩阵，行名是cel id。见(6)。是 'apca' 空间内
        reduction = reduction.norm,  # "pca.norm" 取20个细胞的嵌入均值时，使用的数据名。'pca' 空间内
        dims = 1:ncol(x = embeddings.list.norm[[r]]), #1:PC个数
        verbose = FALSE,
        return.assay = FALSE
      )
    })
    #> cross_impute$pca[[1]][1:2, 1:4] #"pca"表示是 pca 空间内，cross_impute表示使用其他模态最近邻(KNN)挑选的20个细胞
    #         AAACCTGAGACTCGGA-1 AAACCTGAGAGTACAT-1 AAACCTGAGATATACG-1 AAACCTGAGCGAAGGG-1
    #normPC_1        -0.04756215        -0.04552868        -0.07418677          0.5270058
    #normPC_2        -0.75461951        -0.71880930        -0.82398326          0.4956044

    names(x = cross_impute[[r]]) <- cross.modality
    # > str(cross_impute)
    #List of 1
    # $ pca:List of 1 #这里预留空间了，可以支持>=3模态。 在pca空间中，根据 apca KNN挑选的最近邻，求的normPC空间均值
    #  ..$ apca: num [1:30, 1:8036] -0.04756 -0.75462 -0.00974 -0.07171 0.08022 ...
    #  .. ..- attr(*, "dimnames")=List of 2
    #  .. .. ..$ : chr [1:30] "normPC_1" "normPC_2" "normPC_3" "normPC_4" ...
    #  .. .. ..$ : chr [1:8036] "AAACCTGAGACTCGGA-1" "AAACCTGAGAGTACAT-1" "AAACCTGAGATATACG-1" "AAACCTGAGCGAAGGG-1" ...
  }



  #(A14) 模态内预测距离
  within_impute_dist <- lapply(
    X = reduction.list, # list('pca'='pca', 'cpca'='cpca')
    FUN = function(r) { #r='pca' 时
      #1.(原细胞嵌入)x和(最近邻均值估算的细胞嵌入)y的欧氏距离，2.减去 最近一个细胞的距离距离，3.经过 ReLu 负数改为0后，返回
      # 如果是0，表示很近了，如果非0，表示估算位置距离原位置 偏差有点大
      r_dist <- impute_dist( #impute_dist函数见后文(7)。 
        x = query.embeddings.list.norm[[r]], # L2 标准化后的嵌入，默认没有提供，就使用obj自带的细胞嵌入
        #> query.embeddings.list.norm[[r]][1:2, 1:3]
        #                         PC_1       PC_2       PC_3
        #AAACCTGAGACTCGGA-1 0.03717500 -0.7241838 0.16534230
        #AAACCTGAGAGTACAT-1 0.02301086 -0.7453045 0.05282587
        y = t(x = within_impute[[r]]),       # 模态内20个最近邻细胞，求它们在原normPC_空间的位置均值作为估计值
        #> t(within_impute[[r]][1:2, 1:3])
        #                     normPC_1   normPC_2
        #AAACCTGAGACTCGGA-1  0.01897473 -0.8248511
        #AAACCTGAGAGTACAT-1 -0.03725050 -0.7496652
        #AAACCTGAGATATACG-1 -0.06064720 -0.7927484
        nearest.dist = nearest_dist[[r]] #原空间内最近的一个细胞的距离
        # > nearest_dist[[r]] |> head(n=3)
        #AAACCTGAGACTCGGA-1 AAACCTGAGAGTACAT-1 AAACCTGAGATATACG-1 
        #         0.4836937          0.3665269          0.4051788
      )
      return(r_dist) #返回一个 named vector
    }
  )
  #> head(within_impute_dist$pca, n=3)
  #AAACCTGAGACTCGGA-1 AAACCTGAGAGTACAT-1 AAACCTGAGATATACG-1 
  #                 0                  0                  0 


  #(A15) 模态之间预测距离
  cross_impute_dist <- lapply(
    X = reduction.list,
    FUN = function(r) {
      #
      r_dist <- sapply(setdiff(x = reduction.set, y = r), #其他模态都算进去，这里只有一个 'apca'
        FUN = function(r2) { #'apca'
          r2_dist <- impute_dist(
            x = query.embeddings.list.norm[[r]], #原细胞嵌入
            y = t(x = cross_impute[[r]][[r2]]), #使用其他模态的最近邻20个细胞，用它们的原空间位置给出均值作为估计值
            nearest.dist = nearest_dist[[r]] #原空间内最近的一个细胞的距离
          )
          return( r2_dist) #返回一个 named vector
        })
      return(r_dist) #两层结果：输出结果怎么嵌套的？ //见本文3(6)
    }
  )
  #> cross_impute_dist$pca |> head(n=3) #外层是list的键值对，内层是每列一个其他模态
  #                         apca
  #AAACCTGAGACTCGGA-1 0.07634571
  #AAACCTGAGAGTACAT-1 0.00000000
  #AAACCTGAGATATACG-1 0.16821501


  #(A16) 计算核宽度？？
  # calculate kernel width
  #(B1) 实例走这里
  if (snn.far.nn) {
    #(C1) 消息：开始计算核宽度
    if (verbose) {
      message("Calculating kernel bandwidths")
    }
    #(C2) KNN矩阵 获取SNN矩阵
    # 利用细胞之间的最近邻的细胞，求细胞之间的Jacard 相似度 = 交集/并集 * 100%，在[0, 1]之间
    snn.graph.list <- lapply(
      X = sigma.nn.list, # list, 每个模态内20个最近邻的细胞，下标和距离
      FUN = function(nn) { # 遍历每个模态
        snn.matrix <- ComputeSNN( # ComputeSNN 见系列 源码解析21-2(13)，SNN算法见 源码解析21-4(2)
          nn_ranked =  Indices(object = nn)[, 1:s.nn], #给出该模态内20个最近邻cell的下标
          prune = prune.SNN #实例用默认值0
        )
        colnames(x = snn.matrix) <- rownames(x = snn.matrix) <- Cells(x = object)
        return (snn.matrix)
      }
    )
    # > names(snn.graph.list)
    #[1] "pca"  "apca"
    #
    #> snn.graph.list$pca[1:2, 1:3]
    #2 x 3 sparse Matrix of class "dgCMatrix"
    #                   AAACCTGAGACTCGGA-1 AAACCTGAGAGTACAT-1 AAACCTGAGATATACG-1
    #AAACCTGAGACTCGGA-1                  1                  .                  .
    #AAACCTGAGAGTACAT-1                  .                  1                  .

    #(C3) 获取SNN图上距离最远的k.nn个细胞，在原空间L2Norm后的距离的平均值
    farthest_nn_dist <- my.lapply(
      X = 1:length(x = snn.graph.list),
      FUN = function(s) {
        distant_nn <- ComputeSNNwidth(
          snn.graph = snn.graph.list[[s]],
          k.nn = k.nn,
          l2.norm = FALSE,
          embeddings =  embeddings.list.norm[[s]],
          nearest.dist = nearest_dist[[s]]
        )
        return (distant_nn)
      }
    )
    #str(farthest_nn_dist)
    #List of 2
    # $ : num [1:8036] 0.296 0.314 0.34 0.396 0.495 ...
    # $ : num [1:8036] 0.554 0.489 0.506 0.539 0.44 ...

    #(C4) 为list添加名字
    names(x = farthest_nn_dist) <- unlist(x = reduction.list)

    #(C5) 分别乘以1
    modality_sd.list <- lapply(
      X = farthest_nn_dist,
      FUN =  function(sd)  sd * sd.scale
    )
  } else {
  #(B2) 实例不走这里
    if (verbose) {
      message("Calculating sigma by ", sigma.idx, "th neighbor")
    }
    modality_sd.list <- lapply(
      X = reduction.list ,
      FUN =  function(r) {
        rdist <- Distances(object = sigma.nn.list[[r]])[, sigma.idx] - nearest_dist[[r]]
        rdist <- rdist * sd.scale
        return (rdist)
      }
    )
  }


  #(A17) 计算 模态内和之间 的模态核
  # Calculating within and cross modality kernel, and modality weights
  within_impute_kernel <- lapply(
    X = reduction.list,
    FUN = function(r) {
      # 
      exp(-1 * (within_impute_dist[[r]] / modality_sd.list[[r]]) )
    }
  )
  #> within_impute_kernel[[1]] |> head(n=3)
  #AAACCTGAGACTCGGA-1 AAACCTGAGAGTACAT-1 AAACCTGAGATATACG-1 
  #                 1                  1                  1
  # within_impute_kernel[[1]] |> hist(n=100) #[0.5, 1]之间，主要是1。
  # 计算过程: 
  # 1.模态内 最近的20个点的PC均值 / 点到20个最远点距离的均值，
  # 2.分子分母-最近点距离 再取 ReLu
  # y=exp(-x/d)，


  #(A18) 计算 模态之间估算 的模态核
  cross_impute_kernel <- lapply(
    X = reduction.list,
    FUN = function(r) {
      exp(-1 * (cross_impute_dist[[r]] / modality_sd.list[[r]]) )
    }
  )
  #> cross_impute_kernel[[1]] |> head(n=2)
  #                        apca
  #AAACCTGAGACTCGGA-1 0.7729464
  #AAACCTGAGAGTACAT-1 1.0000000
  # cross_impute_kernel[[1]][, 'apca'] |> hist(n=100) #分布到[0,1]之间，0.5最多
  # 计算过程: 
  # 1.模态间 最近的20个点的PC均值 / 点到20个最远点距离的均值，
  # 2.分子分母-最近点距离 再取 ReLu
  # y=exp(-x)，


  #(A19) 记录参数，写入日志？
  params <- list(
    "reduction.list" = reduction.list,
    "dims.list" = dims.list,
    "l2.norm" = l2.norm,
    "k.nn" = k.nn,
    "sigma.idx" = sigma.idx,
    "snn.far.nn" = snn.far.nn ,
    "sigma.list" = modality_sd.list,
    "nearest.dist" = nearest_dist
  )


  #(A20) 模态打分： 就是求比值， 模态内核 / 模态之间核
  modality_score <-  lapply(
    X = reduction.list,
    FUN = function(r) {
      score.r <- sapply(
        X = setdiff(x = reduction.set, y = r),
        FUN = function(r2) {
          # 组内 / (组间 + 1e-4)
          score <- within_impute_kernel[[r]] / (cross_impute_kernel[[r]][, r2] + cross.contant.list[[r]])
          # 低于min的设为min，超过max的设为max
          score <- MinMax(data = score, min = 0, max = 200) #MinMax 见本系列 源码解析19-2(9)
          return(score)
        }
      )
      return(score.r)
    }
  )
  # > dim(modality_score[[1]])  #[1] 8036    1
  #> modality_score[[1]][1:3, 1, drop=F]
  #                       apca
  #AAACCTGAGACTCGGA-1 1.293583
  #AAACCTGAGAGTACAT-1 0.999900
  #AAACCTGAGATATACG-1 1.640346
  # modality_score[[1]][,1] |> hist(n=100) #[0, 50] 主峰是[0,10]


  #(A21) 实例是F，跳过
  if (smooth) {
    # But: 还是看看怎么进行平滑化的: 用最近20个细胞的 模态平均值
    modality_score <- lapply(
      X = reduction.list,
      FUN = function(r) {
        apply(
          X = Indices(object = nn.list[[r]]), #每行是1个细胞的20个最近邻cell的下标
          #> str(nn.list[['pca']]@nn.idx)
          # num [1:8036, 1:20] 1 2 3 4 5 6 7 8 9 10 ...
          #> nn.list[['pca']]@nn.idx[1:2, 1:4]
          #     [,1] [,2] [,3] [,4]
          #[1,]    1 4183 6473 5855
          #[2,]    2 3214  141 4925
          MARGIN = 1,
          FUN = function(nn)  mean(x = modality_score[[r]][nn[-1]]) #这些cell除了自己，其余求平均值
          # 如果是多于3个组学，这里可能有问题: 无法保持原来的维度。
          # > modality_score[['pca']] [1:2, 1, drop=F]
          #                       apca
          #AAACCTGAGACTCGGA-1 1.293583
          #AAACCTGAGAGTACAT-1 0.999900
          # > modality_score[['pca']] [1:2, 1, drop=F] |> class()
          #[1] "matrix" "array"
          #
          # 例: 只取出第一列
          # > t1=iris[1:4, 1:3] |> as.matrix(); t1
          #  Sepal.Length Sepal.Width Petal.Length
          #1          5.1         3.5          1.4
          #2          4.9         3.0          1.4
          #3          4.7         3.2          1.3
          #4          4.6         3.1          1.5
          #> class(t1)
          #[1] "matrix" "array" 
          #> t1[c(1,3,4)]
          #[1] 5.1 4.7 4.6
          #
          # 应该加逗号，和 drop=F：经测试，这样也不行，同样无法保持 modality_score 的维度经过smooth不变...//todo
          # > t1[c(1,3,4), , drop=F]
          #  Sepal.Length Sepal.Width Petal.Length
          #1          5.1         3.5          1.4
          #3          4.7         3.2          1.3
          #4          4.6         3.1          1.5
        )
      }
    )
  }
  #> modality_score3$pca[1:3]
  #AAACCTGAGACTCGGA-1 AAACCTGAGAGTACAT-1 AAACCTGAGATATACG-1 
  #          1.565276           1.235708           1.584979 


  #(A22) softmax的分母: e(x1)/[e(x1) + e(x2) + ...]
  # 所有模态打分，求exp()后，按行求和
  all_modality_score <- rowSums(x = exp(x = Reduce(f = cbind, x = modality_score)))
  #> all_modality_score[1:2]
  #AAACCTGAGACTCGGA-1 AAACCTGAGAGTACAT-1
  #          3.645828           2.718010

  #(A23) softmax得到模态权重：exp(模态1打分) / 累加[exp(模态1打分, exp(模态2打分), ...)]
  # 卷积神经网络常用，不多说。
  modality.weight <- lapply(
    X = modality_score,
    FUN = function(score_m) {
      rowSums(x = exp(x = score_m))/all_modality_score
    }
  )
  #> modality.weight$pca[1:2]
  #AAACCTGAGACTCGGA-1 AAACCTGAGAGTACAT-1 
  #                 1                  1
  # 使用全部数据时，结果如下：
  #> modality.weight_2@modality.weight.list$pca[1:2]
  #AAACCTGAGACTCGGA-1 AAACCTGAGAGTACAT-1 
  #         0.5325040          0.4924699 


  #(A24) 记录重要的中间结果，用于构建返回值
  score.mat <- list(
    within_impute_dist = within_impute_dist, #模态内KNN最近邻20个cell，在原normPC的坐标均值
    cross_impute_dist = cross_impute_dist, #其他模态KNN给出的最近邻20个cell，在原normPC的坐标均值
    within_impute_kernel = within_impute_kernel, #模态内SNN给出的最远的20个非零cell，在原normPC距离该细胞的距离，均值
    cross_impute_kernel = cross_impute_kernel, #其他模态SNN给出的最远的20个非零cell，在原normPC距离该细胞的距离，均值
    modality_score = modality_score #模态打分，exp前
  )


  #(A25) 日志_Seurat命令
  # unlist the input parameters
  command <- LogSeuratCommand(object = object, return.command = TRUE)
  command@params <- lapply(X = command@params, FUN = function(l) unlist(x = l))


  #(A26) 模态用到的实验名字
  modality.assay <- sapply(
    X = reduction.list ,
    FUN = function (r) slot(object[[r]], name = "assay.used")
  )
  #> object[["pca"]]@assay.used
  #[1] "RNA"


  #(A27) 返回值
  modality.weights.obj <- new(
    Class = "ModalityWeights", #模态权重 对象
    modality.weight.list = modality.weight, #权重，最希望的结果
    modality.assay = modality.assay, #模态用的实验名
    params = params, #参数
    score.matrix = score.mat, #打分矩阵
    command = command
  )

  return(modality.weights.obj)
}








(2B) FindModalityWeights() 第二遍阅读：代码块阅读

目的：由结果反推过程。修正错误。

```
> names(bm@assays)
[1] "RNA" "ADT"

> names(bm@reductions)
[1] "pca"  "apca"

> bm@reductions$pca@assay.used
[1] "RNA"
> bm@reductions$apca@assay.used
[1] "ADT"

> bm
An object of class Seurat 
36656 features across 8036 samples within 2 assays 
Active assay: RNA (36601 features, 2000 variable features)
 1 other assay present: ADT
 2 dimensional reductions calculated: pca, apca


主函数：
> modality.weight_2 <- Seurat:::FindModalityWeights(
  object = bm,
  reduction.list = list("pca", "apca"),
  dims.list = list(1:30, 1:20),
  k.nn = 20,
  sd.scale = 1,
  l2.norm = T,
  cross.contant.list = NULL,
  smooth = F,
  verbose = T
)
Finding 20 nearest neighbors for each modality.
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=01m 16s
Calculating kernel bandwidths
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=29s

> str(modality.weight_2) #输出太多，不列举了

> slotNames(modality.weight_2)
[1] "modality.weight.list" "modality.assay"       "params"               "score.matrix"         "command"

其中
"modality.weight.list"：模态权重，期望的输出。 2个list，pca和apca
"modality.assay": 2个字符串，"RNA", "ADT"
"params": 输入参数
"score.matrix": 重要的中间过程，我们接着代码块核对的重点
"command": 命令记录


> modality.weight_2@score.matrix |> names()
[1] "within_impute_dist"   "cross_impute_dist"    "within_impute_kernel" "cross_impute_kernel"  "modality_score"
```

1. impute细胞：按照模态内和外的KNN给出20个最近邻细胞，在原normPC中计算各维度平均值。
2. 预测距离：impute的细胞 与 原细胞，在NormPC中的欧氏距离，减去其与原细胞最近一个细胞的距离，取ReLu。
3. 核宽度 kernel bandwidths: 取SNN图中与原细胞最远的20个细胞，与原细胞在NormPC中的欧氏距离，减去最近一个细胞距离，取Relu
4. 模态内/外预测核(impute kernel)： epx(- 模态内预测距离 / 模态宽度 )
5. 模态打分(modality score)： 模态内预测核 / (模态间预测核+1e-4)
6. softmax得到模态权重(modality weight): exp(打分1) / 累加(exp(打分1), exp(打分2), ...)



=> 该问题依旧：能处理>=3个模态吗？我倾向于认为可以。
1. 代码实现倾向于支持多个模态间的列。
2. 最后是softmax 把打分转为比例，且比例之和为1，这一步也支持多分类。






(2C) FindModalityWeights() 第三遍阅读：结合外部资料，思考WNN全过程
# https://www.cell.com/cell/fulltext/S0092-8674(21)00583-3
# https://mp.weixin.qq.com/s?__biz=MjM5NTk0Mzg2Nw==&mid=2247489705&idx=1&sn=a726e941173f10eb96e82168c1d7a817
# https://blog.csdn.net/weixin_53637133/article/details/138020958

每个模态内部按Seurat标准流程做到PCA，然后用PCA获得的KNN计算其模态的权重。

WNN的构建过程分为四步：
（1）Constructing independent knearest neighbor (KNN) graphs for both modalities。
（2）Performing within and across-modality prediction
* within-modality predictions：对于细胞中每个RNA而言，根据RNA数据中K个最临近细胞中的RNA值取均值得到预测的RNA值；对于细胞中每个蛋白而言，根据蛋白数据中K个最临近细胞中的蛋白值取均值得到预测的蛋白值。
* cross-modality predictions：对于细胞中每个RNA而言，根据蛋白数据中K个最临近细胞中的RNA值取均值得到预测的RNA值；对于细胞中每个蛋白而言，根据RNA数据中K个最临近细胞中的蛋白值取均值得到预测的蛋白值。
（3）Calculating cell-specific modality weights: 根据预测值和实际观测值，利用UMAP算法中的指数核函数计算亲和度，并使用softmax转换得到针对蛋白和RNA的权重值。
（4）Calculating a WNN graph: 根据权重值，重建考虑权重后的KNN图，即WNN图


i) bandwidth的定义主要依据该细胞与邻居之间欧氏距离的平均值。

ii) 每次求距离，都要减去原细胞最近一个细胞的距离，为什么？
有论文（Umap: Uniform manifold approximation and projection for dimension reduction）证明：
一个细胞和它最临近细胞之间的距离，会使得计算过程中分子和分母出现膨胀问题，所以我们需要在计算过程中对分子和分母都减去一个细胞和它最临近细胞之间的距离














(3) L2Norm() 对矩阵进行L2标准化：每行除以该行的 平方和的开方
$ find . | grep -v "dustbin/" | grep -v "notes/" | xargs grep -in "L2Norm" --color=auto 2>/dev/null
./seurat-4.4.0/R/dimensional_reduction.R:2256:L2Norm <- function(vec) { #这个是对向量的
./seurat-4.4.0/R/utilities.R:1720:L2Norm <- function(mat, MARGIN = 1){ #像是这个，对矩阵的

问题：共2个同名函数，原来的包怎么保证该使用哪个呢？这里没有OOP机制。
默认暴漏的是矩阵版 `> Seurat:::L2Norm`。可能是因为它是后来定义的，覆盖了前者？ //todo


# Calculates the l2-norm of a vector 对向量
# 这个实现很奇怪，只是平方和的开方，就是欧氏距离。没有做除法，这叫标准化吗？//todo
#
# Modified from PMA package
# @references Witten, Tibshirani, and Hastie, Biostatistics 2009
# @references \url{https://github.com/cran/PMA/blob/master/R/PMD.R}
#
# @param vec numeric vector
#
# @return returns the l2-norm.
#
L2Norm <- function(vec) {
  a <- sqrt(x = sum(vec ^ 2)) #平方和再开方
  if (a == 0) { #如果等于0，则改为0.5
    a <- .05
  }
  return(a)
}




# L2 normalize the columns (or rows) of a given matrix | 对矩阵的列（或行）做L2标准化
# @param mat Matrix to cosine normalize | 将要做cosine标准化的矩阵
# @param MARGIN Perform normalization over rows (1) or columns (2) | 默认对行做标准化
#
# @return returns l2-normalized matrix
#
L2Norm <- function(mat, MARGIN = 1){
  normalized <- Sweep(
    x = mat,
    MARGIN = MARGIN, #默认，按行。除以每一行的
    STATS = apply( # 行：平方和的开方
      X = mat,
      MARGIN = MARGIN,
      FUN = function(x){
        sqrt(x = sum(x ^ 2))
      }
    ),
    FUN = "/"
  )
  # 无穷大都改写为0： 除数是0导致无穷大，这一行全是0导致的
  normalized[!is.finite(x = normalized)] <- 0
  return(normalized)
}


# 手工实现
tmp1=sweep(embeddings.list$pca,
           MARGIN = 1, 
           STATS = apply(embeddings.list$pca, 1, function(x){
            sqrt( sum(x**2) )
          }),
          FUN = "/"
)
# 调用函数
tmp2 = Seurat:::L2Norm(embeddings.list$pca)

all(abs(tmp1-tmp2)<1e-10) #T
tmp1[1:3, 1:4]
tmp2[1:3, 1:4]
#                          PC_1       PC_2        PC_3        PC_4
#AAACCTGAGACTCGGA-1  0.03717500 -0.7241838  0.16534230  0.06183454
#AAACCTGAGAGTACAT-1  0.02301086 -0.7453045  0.05282587  0.04039307
#AAACCTGAGATATACG-1 -0.05451243 -0.7092048 -0.01946849 -0.02317811

问题：输入是行为cell id，列为 PC_1 的cell embedding 矩阵，为什么对行做标准化？ 
我的解释：保证每一行的平方和都是1，每个点都在距离原点单位1半径的球面上。











(4) Distances() 从 "Neighbor"对象返回带cell id的距离矩阵

r="pca"
t1=Distances(object = nn.list[[r]])

> class(t1)
[1] "matrix" "array" 

输入：
> nn.list[['pca']]@nn.dist[1:3, 1:4]
     [,1]      [,2]      [,3]      [,4]
[1,]    0 0.4836937 0.4968935 0.4984034
[2,]    0 0.3665269 0.3854056 0.3913714
[3,]    0 0.4051788 0.4284064 0.4461618

> t1[1:3, 1:4]
                   [,1]      [,2]      [,3]      [,4]
AAACCTGAGACTCGGA-1    0 0.4836937 0.4968935 0.4984034
AAACCTGAGAGTACAT-1    0 0.3665269 0.3854056 0.3913714
AAACCTGAGATATACG-1    0 0.4051788 0.4284064 0.4461618          

这个函数仅仅是添加cell id? 看了源码，确实是返回带cell id的距离矩阵。



# 第一列是自己，第二列是最近的邻居
nearest_dist <- lapply(X = reduction.list, FUN = function(r) Distances(object = nn.list[[r]])[, 2])
> str(nearest_dist)
#List of 2
# $ pca : Named num [1:8036] 0.484 0.367 0.405 0.276 0.625 ...
#  ..- attr(*, "names")= chr [1:8036] "AAACCTGAGACTCGGA-1" "AAACCTGAGAGTACAT-1" "AAACCTGAGATATACG-1" "AAACCTGAGCGAAGGG-1" ...
# $ apca: Named num [1:8036] 0.452 0.452 0.555 0.561 0.827 ...
#  ..- attr(*, "names")= chr [1:8036] "AAACCTGAGACTCGGA-1" "AAACCTGAGAGTACAT-1" "AAACCTGAGATATACG-1" "AAACCTGAGCGAAGGG-1" ...




$ find . | grep -v "dustbin/" | grep -v "notes/" | xargs grep -n "Distances" --color=auto 2>/dev/null
./seurat-object-4.1.4/NAMESPACE:85:S3method(Distances,Neighbor) #这是一个S3类

#./seurat-object-4.1.4/R/generics.R:458:Distances <- function(object, ...) {
#./seurat-object-4.1.4/R/generics.R:459:  UseMethod(generic = 'Distances', object = object)

#' Get the Neighbor nearest neighbors distance matrix
#'
#' @param object An object
#' @param ... Arguments passed to other methods
#'
#' @return The distance matrix
#'
#' @rdname Distances
#' @export Distances
#'
#' @concept data-access
#'
Distances <- function(object, ...) {
  UseMethod(generic = 'Distances', object = object)
}


#./seurat-object-4.1.4/R/neighbor.R:74:Distances.Neighbor <- function(object, ...) {

#' @rdname Distances
#' @export
#' @method Distances Neighbor
#'
Distances.Neighbor <- function(object, ...) {
  object <- UpdateSlots(object = object)
  # 获取距离
  distances <- slot(object = object, name = "nn.dist")
  # 添加行名：细胞名
  rownames(x = distances) <- slot(object = object, name = "cell.names")
  return(distances)
}










(5) PredictAssay() 最近邻的19个细胞（去掉自己）的均值，是其在 normPC_x 空间的预测值

$ find . | grep -v "dustbin/" | grep -v "notes/" | xargs grep -n "PredictAssay" --color=auto 2>/dev/null
./seurat-4.4.0/R/clustering.R:210:PredictAssay <- function(
./seurat-4.4.0/R/clustering.R:1125:    within_impute[[r]] <- PredictAssay(
./seurat-4.4.0/R/clustering.R:1135:      PredictAssay(


#' Predict value from nearest neighbors | 从最近邻预测值
#'
#' This function will predict expression or cell embeddings from its k nearest neighbors index. | 该函数将从其k最近邻预测其表达量或者细胞嵌入
#' For each cell, it will average its k neighbors value to get its new imputed value. | 对每个细胞，取k最近邻的均值作为新的预测值
#' It can average expression value in assays and cell embeddings from dimensional reductions. | 可从降维中的assays和细胞嵌入获取表达平均值
#'
#' @param object The object used to calculate knn | 用于计算knn的对象
#' @param nn.idx k near neighbour indices. A cells x k matrix. | k最近邻的索引：cell * k 矩阵
#' @param assay Assay used for prediction | 用于预测的assay
#' @param reduction Cell embedding of the reduction used for prediction | 用于预测的降维后的细胞嵌入
#' @param dims Number of dimensions of cell embedding | 细胞嵌入的维度数
#' @param return.assay Return an assay or a predicted matrix | 返回assay还是预测的矩阵
#' @param slot slot used for prediction | 用于预计的slot，默认是 scale.data
#' @param features features used for prediction | 用于预测的基因
#' @param mean.function the function used to calculate row mean | 用于计算行平均值的函数
#' @param seed Sets the random seed to check if the nearest neighbor is query cell | 随机数种子，检查最近邻是否是查询的细胞
#' @param verbose Print progress | 是否打印进度条
#'
#' @return return an assay containing predicted expression value in the data slot | 返回assay，data slot 包含预测表达值
#' @concept integration
#' @export
#'
PredictAssay <- function(
  object,
  nn.idx,
  assay,
  reduction = NULL,
  dims = NULL,
  return.assay = TRUE,
  slot = "scale.data",
  features = NULL,
  mean.function = rowMeans,
  seed = 4273,
  verbose = TRUE #函数内调用时，设置为F了
){
  
  #(A1) 如果 行均值 函数不是函数类型，则报错
  if (!inherits(x = mean.function, what = 'function')) {
    stop("'mean.function' must be a function")
  }

  #(A2) 如果 reduction 为空：本实例不为空，为 "pca.norm"
  if (is.null(x = reduction)) {
  	#(B1) 实例跳过本if
    reference.data <- GetAssayData(
      object = object,
      assay = assay,
      slot = slot
    )
    features <- features %||% VariableFeatures(object = object[[assay]])
    if (length(x = features) == 0) {
      features <- rownames(x = reference.data)
      if (verbose) {
        message("VariableFeatures are empty in the ", assay,
                " assay, features in the ", slot, " slot will be used" )
      }
    }
    reference.data <- reference.data[features, , drop = FALSE]
  } else {
  	#(B2) 实例走else
  	# 如果dims为空，则报错。实例不空
    if (is.null(x = dims)) {
      stop("dims is empty")
    }
    # 参考数据：获取 "pca.norm" 的 1:xx 列； 矩阵转置
    reference.data <- t(x = Embeddings(object = object, reduction = reduction)[, dims])
    #         AAACCTGAGACTCGGA-1 AAACCTGAGAGTACAT-1
	#normPC_1          0.0371750         0.02301086
	#normPC_2         -0.7241838        -0.74530448
  }

  #(A3) 设置随机数：仅仅用于检测第一列是否是细胞自身的下标
  set.seed(seed = seed)
  #(A4) 抽样。 样本是1:行细胞个数，抽样大小是 min(50, 细胞数)
  nn.check <- sample(x = 1:nrow(x = nn.idx), size = min(50, nrow(x = nn.idx)))
  #(A5) 如果第一列是细胞自己，去掉第一列
  if (all(nn.idx[nn.check, 1] == nn.check)) {
    if(verbose){
      message("The nearest neighbor is the query cell itself, and it will not be used for prediction")
    }
    nn.idx <- nn.idx[,-1] #行：该细胞的20-1=19个最近邻的cell索引下标
    #                   [,1] [,2] [,3] [,4]
	#AAACCTGAGACTCGGA-1 4183 6473 5855 4815
	#AAACCTGAGAGTACAT-1 3214  141 4925 6424
  }

  #(A6) 按行(cell id)求：每个细胞的 19个最近邻细胞的 normPC_x 的行平均值
  # 仔细品味该apply细节，其返回值也是本函数的 返回值。
  # 输入x是行 nn.idx 的行，行名是 cell id，列为其最近邻cell的index。
  # 按行，取出该cell的最近邻细胞下标。
  # ref 结构：行为 normPC_d, 列为 cell id。
  # 按列取ref中这些最近邻细胞列index，求行平均，获得每个 normPC_ 的均值。
  # 就是这个cell id 根据 组学最近邻(nn.idex)，在原 ref空间 的normPC_中的预测位置
  predicted <- apply( 
    X = nn.idx,
    MARGIN = 1,
    FUN = function(x) mean.function(reference.data[, x] )
  )
  #设置行名：cell id
  colnames(x = predicted) <- Cells(x = object)

  #(A7) 调用时设置为F
  if (return.assay) {
  	#(B1)实例不走这
  	# 创建一个assay对象，其 data slot 是预测结果
    predicted.assay <- CreateAssayObject(data = predicted, check.matrix = FALSE)
    return (predicted.assay)
  } else {
  	#(B2)实例走这，直接返回(A6)矩阵
    return (predicted)
    #> predicted[1:2, 1:2]
	#         AAACCTGAGACTCGGA-1 AAACCTGAGAGTACAT-1
	#normPC_1         0.01897473         -0.0372505
	#normPC_2        -0.82485107         -0.7496652
	#> dim(predicted)
	#[1]   30 8036
  }

}








(6) Indices() 返回最近邻索引 cell*k 矩阵，行名是cel id，值是最近邻cell下标
#./seurat-object-4.1.4/R/generics.R:704:Indices <- function(object, ...) {

#' Get Neighbor nearest neighbor index matrices #泛型函数
#'
#' @param object An object
#' @param ... Arguments passed to other methods;
#'
#' @return A matrix with the nearest neighbor indices
#'
#' @rdname Indices
#' @export Indices
#'
#' @concept data-access
#'
Indices <- function(object, ...) {
  UseMethod(generic = "Indices", object = object)
}



#./seurat-object-4.1.4/R/neighbor.R:110:Indices.Neighbor <- function(object, ...) {


#' @rdname Indices
#' @export
#' @method Indices Neighbor
#'
Indices.Neighbor <- function(object, ...) {
  object <- UpdateSlots(object = object)
  indices <- slot(object = object, name = "nn.idx") #返回其 nn.idex，就是最近邻索引 cell*k 矩阵
  rownames(x = indices) <- slot(object = object, name = "cell.names") #行名是cell id
  return(indices)
}






(7) impute_dist() 计算 x 和 y 的欧式距离 - 最近一个细胞的距离，取ReLu
$ find . | grep -v "dustbin/" | grep -v "notes/" | xargs grep -n "impute_dist" --color=auto 2>/dev/null

# ./seurat-4.4.0/R/utilities.R:2040:impute_dist <- function(x, y, nearest.dist) {

# Calculate euclidean distance the x and y, | 计算x和y的欧氏距离
# and subtract the nearest neighbors of x distance to keep local connectivity | 减去最近一个邻居的距离，保持局部可连接性
# It is used in FindModalityWeights to calculate the with and cross modality distance | FindModalityWeights() 用于计算模态内和模态之间的距离
impute_dist <- function(x, y, nearest.dist) {
  #(1) 计算每个细胞在两个模态中距离：
  # 行：（x和y之间的欧氏距离）差的平方和的开方 - 距离最近的一个细胞的距离
  dist <- sqrt(x = rowSums(x = (x - y)**2)) - nearest.dist
  #(2) 只保留非负项
  dist <- ReLu(x = dist)
  return(dist)
}

# 








(8) ReLu() 负数都改为0
# Rectified linear units function. Calculate positive part of its argument
# The input can be a vector and a matrix
ReLu <- function(x) {
  x[x < 0] <- 0
  return(x)
}

# 








(9) ComputeSNNwidth()：每个细胞与SNN中的k.nn个最远邻居的平均距离

$ find . | grep -v "dustbin/" | grep -v "notes/" | xargs grep -n "ComputeSNNwidth" --color=auto 2>/dev/null
./seurat-4.4.0/R/clustering.R:966:ComputeSNNwidth <- function(
./seurat-4.4.0/R/clustering.R:1192:        distant_nn <- ComputeSNNwidth(


# Calculate mean distance of the farthest neighbors from SNN graph | 计算SNN图中最远邻居的平均距离
#
# This function will compute the average distance of the farthest k.nn neighbors with the lowest nonzero SNN edge weight. 
# 计算 非零SNN边权重 的 最远 k.nn 邻居平均距离。
# First, for each cell it finds the k.nn neighbors with the smallest edge weight. 
# 每个细胞，在20个最近邻中，找最小边权重
# If there are multiple cells with the same edge weight at the k.nn-th index, consider all of those cells in the next step.
# 如果同一个边权重上有多个细胞下标，它们都进入下一步。
#
# Next, it computes the euclidean distance to all k.nn cells in the space defined by the embeddings matrix 
# 接着，在细胞嵌入的空间中，计算其与k.nn个细胞的欧式距离。
# 
# and returns the average distance to the farthest k.nn cells.
# 然后，返回到最远k.nn细胞的平均距离
#
# @param snn.graph An SNN graph | 一个SNN图
# @param embeddings The cell embeddings used to calculate neighbor distances | 用于计算相邻距离的细胞嵌入
# @param k.nn The number of neighbors to calculate | 用于计算的邻居个数
# @param l2.norm Perform L2 normalization on the cell embeddings | 是否在细胞嵌入上进行L2标准化：默认T
# @param nearest.dist The vector of distance to the nearest neighbors to subtract off from distance calculations
# 到最近的细胞的距离，从欧氏距离中减去。
#
#
ComputeSNNwidth <- function(
  snn.graph,
  embeddings,
  k.nn,
  l2.norm = TRUE,
  nearest.dist = NULL
) {

  #(A1) 实例不做 l2.norm=F
  if (l2.norm) {
    embeddings <- L2Norm(mat = embeddings)
  }

  #(A2) 实例不是空
  nearest.dist <- nearest.dist %||% rep(x = 0, times = ncol(x = snn.graph))

  #(A3) list长度不同，则报错
  if (length(x = nearest.dist) != ncol(x = snn.graph)) {
    stop("Please provide a vector for nearest.dist that has as many elements as",
         " there are columns in the snn.graph (", ncol(x = snn.graph), ").")
  }

  #(A4) 就这一个干活的：c++函数，见(10)
  snn.width <- SNN_SmallestNonzero_Dist(
    snn = snn.graph, #SNN图
    mat = embeddings, #L2标准化的细胞嵌入
    n = k.nn, #20
    nearest_dist = nearest.dist #最近一个细胞的距离
  )
  return (snn.width)
}









(10) SNN_SmallestNonzero_Dist() 是一个c++函数：每个细胞与SNN上最远的点的欧氏距离的平均值

* 每个细胞的SNN非0且相似性最小的n个细胞
* 与该细胞在原空间L2标准化的距离-最近一个细胞的距离
* 去掉负数后，求均值

$ find . | grep -v "dustbin/" | grep -v "notes/" | xargs grep -n "SNN_SmallestNonzero_Dist" --color=auto 2>/dev/null
./seurat-4.4.0/src/snn.cpp:82:std::vector<double> SNN_SmallestNonzero_Dist(


// [[Rcpp::export]]
std::vector<double> SNN_SmallestNonzero_Dist(
    Eigen::SparseMatrix<double> snn, //稀疏矩阵 snn 相似性图，对称矩阵
    Eigen::MatrixXd mat,            //矩阵 mat：细胞嵌入图，L2标准化后的
    int n,      //实例为 20
    std::vector<double> nearest_dist //最近一个细胞的距离 浮点数向量
) {
  # 注： c++的注释//写法在当前文件中不会被识别为注释，同一个文件中，当前能识别的注释只有#
  #(A1) 声明变量 results，作为返回值
  std::vector<double> results;
  
  #(A2) 对矩阵外循环(是 列数) 
  # snn 行和列都是cell id, 值是相互之间的Jacard 相似度。 SNN图是对称矩阵，行和列都一样。
  for (int i=0; i < snn.outerSize(); ++i){
    # > snn[1:2, 1:2]
    #2 x 2 sparse Matrix of class "dgCMatrix"
    #                   AAACCTGAGACTCGGA-1 AAACCTGAGAGTACAT-1
    #AAACCTGAGACTCGGA-1                  1                  .
    #AAACCTGAGAGTACAT-1                  .                  1

    #(B1) 声明向量，用于保存非零snn的元素，和下标
    // create vectors to store the nonzero snn elements and their indices
    std::vector<double> nonzero;
    #(B2) 见上：用于保存 非零 下标
    std::vector<size_t> nonzero_idx;

    #(B3) 对该列的 非零元素 遍历。保存值和行号
    for (Eigen::SparseMatrix<double>::InnerIterator it(snn, i); it; ++it) {
      nonzero.push_back(it.value());
      nonzero_idx.push_back(it.row());
    }

    #(B4) 对值进行排序：测试是 升序，从小到大 // sort_indexes() 见下文(11)
    std::vector<size_t> nonzero_order = sort_indexes(nonzero); # 这个排序有意义吗？接着是遍历所有非零值，求距离的平均

    #(B5) 最近邻个数
    int n_i = n;

    #(B6) n_i 等于实际非0值个数
    if (n_i > nonzero_order.size()) n_i = nonzero_order.size();

    #(B7) 向量 dists，记录距离
    std::vector<double> dists;

    #(B8) 遍历非零值 向量: （外循环一列中）从snn相似性最低 到 最高
    for (int j = 0; j < nonzero_order.size(); ++j) {
      // compute euclidean distances to cells with small edge weights
      // if multiple entries have same value as nth element, calc dist to all
      # 获取一个值 对应的 cell 下标
      size_t cell = nonzero_idx[nonzero_order[j]];

      #(C1) 如果dists长度小于 n， 或者 非零值 等于 j 和 n-1 是同一个（最后一个元素了）
      if(dists.size() < n_i  || nonzero[nonzero_order[j]] == nonzero[nonzero_order[n_i-1]]) {

        #(D1) 细胞嵌入图 mat: 外循环i(列名)细胞 和 内循环(本列非零值)cell细胞 差，的 norm() 方法
        # mat[1:2, 1:4] # mat 就是 embeddings.list.norm[['pca']]，是L2标准化后的
        #                         PC_1       PC_2       PC_3       PC_4
        #AAACCTGAGACTCGGA-1 0.03717500 -0.7241838 0.16534230 0.06183454
        #AAACCTGAGAGTACAT-1 0.02301086 -0.7453045 0.05282587 0.04039307
        double res = (mat.row(cell) - mat.row(i)).norm(); #.norm()是求L2范数（欧氏距离）：差的平方和，再开方

        #(D2) 如果原空间中最近一个细胞距离 >0
        if (nearest_dist[i] > 0) {
          # 则 刚计算的L2范数 - 原最近距离
          res = res - nearest_dist[i];
          # 取ReLu：负数改为0
          if (res < 0) res = 0;
        }
        # 作为相对欧氏距离，压入dists向量末尾 
        dists.push_back(res);


      #(C2) 否则，内循环结束。什么情况下发生？
      # 我没想到
      } else {
        break;
      }
    }

    #(B9) 声明变量：平均距离
    double avg_dist;

    #(B10) 距离累加成一个数，除以个数，得距离均值
    if (dists.size() > n_i) {
      std::sort(dists.rbegin(), dists.rend());
      avg_dist = std::accumulate(dists.begin(), dists.begin() + n_i, 0.0) / n_i;
    } else {
      avg_dist = std::accumulate(dists.begin(), dists.end(), 0.0) / dists.size();
    }

    #(B11) 获得该细胞的距离均值
    results.push_back(avg_dist);
  }

  #(A3)
  return results;
}










(11) sort_indexes() C++函数：和R的order()效果一致，只是R是1-based下标

$ find . | grep -v "dustbin/" | grep -v "notes/" | xargs grep -n "sort_indexes" --color=auto 2>/dev/null
./seurat-4.4.0/src/data_manipulation.cpp:416:std::vector<size_t> sort_indexes(const std::vector<S> &v) {
./seurat-4.4.0/src/integration.cpp:99:std::vector<size_t> sort_indexes(const std::vector<S> &v) {
./seurat-4.4.0/src/snn.cpp:72:std::vector<size_t> sort_indexes(const std::vector<S> &v) {
./seurat-object-4.1.4/src/data_manipulation.cpp:15:std::vector<size_t> sort_indexes(const std::vector<S> &v) {
这个函数一模一样的找到4份。

template <typename S>
std::vector<size_t> sort_indexes(const std::vector<S> &v) {
  // initialize original index locations
  std::vector<size_t> idx(v.size()); //声明空数组，等长度
  std::iota(idx.begin(), idx.end(), 0); //从0开始递增填充：0,1,2,...
  std::stable_sort(idx.begin(), idx.end(), //
                   [&v](size_t i1, size_t i2) {return v[i1] < v[i2];});
  return idx;
}



测试：
$ cat test_mat6.cpp
#include <Rcpp.h>
using namespace Rcpp;

template <typename S>
std::vector<size_t> sort_indexes(const std::vector<S> &v) {
  // initialize original index locations
  std::vector<size_t> idx(v.size()); //声明空数组，等长度
  std::iota(idx.begin(), idx.end(), 0); //从0开始递增填充：0,1,2,...
  std::stable_sort(idx.begin(), idx.end(), //
                   [&v](size_t i1, size_t i2) {return v[i1] < v[i2];});
  return idx;
}

// [[Rcpp::export]]
void fn6() {
  // accumulate?
  std::vector<double> nonzero={2.2, 1.1, 3.3, 1.1, 4.4};
  Rcpp::Rcout << "nonzero.size():" << nonzero.size() <<  std::endl;
  
  
  std::vector<size_t> nonzero_order = sort_indexes(nonzero);
  for(int i=0; i<nonzero_order.size(); i++){
    size_t idx = nonzero_order[i];
    Rcpp::Rcout << i << " " << idx <<" " << nonzero[idx] << std::endl;
  }
}


> Rcpp::sourceCpp("test_mat6.cpp", cacheDir="./tmp/", verbose=T, rebuild=T)

> fn6()
nonzero.size():5
0 1 1.1
1 3 1.1
2 0 2.2
3 2 3.3
4 4 4.4

结论： sort_indexes() 返回值 是输出第二列，可见是原数组的下标。
  其中两个1.1出现2次，按出现顺序排列。

> a1=c(2.2, 1.1, 3.3, 1.1, 4.4)
> order(a1)
[1] 2 4 1 3 5
> order(a1)-1 #和R的order()函数的效果一致，只是R是1-based下标。
[1] 1 3 0 2 4










(12) MultiModalNN() 找多模态最近邻
$ find . | grep -v "dustbin/" | grep -v "notes/" | xargs grep -n "MultiModalNN" --color=auto 2>/dev/null
./seurat-4.4.0/R/clustering.R:1379:MultiModalNN <- function(


# Find multimodal neighbors
#
# @param object The object used to calculate knn | 计算knn的obj
# @param query The query object when query and reference are different | 查询对象
# @param modality.weight A \code{\link{ModalityWeights}} object generated by
# \code{\link{FindModalityWeights}} | FindModalityWeights() 返回的模态权重对象
# @param modality.weight.list A list of modality weight value | 模态权重的值
# @param k.nn Number of nearest multimodal neighbors to compute | 要计算的多模态最近邻个数
# @param reduction.list A list of reduction name | 降维名字
# @param dims.list A list of dimensions used for the reduction |  降维的维度
# @param knn.range The number of approximate neighbors to compute | 要计算的最近邻个数
# @param kernel.power The power for the exponential kernel | 指数核函数的幂
# @param nearest.dist The list of distance to the nearest neighbors | 到最近邻的距离
# @param sigma.list The list of kernel width | 核函数宽度
# @param l2.norm Perform L2 normalization on the cell embeddings after
# dimensional reduction | 是否对细胞嵌入做L2标准化
# @param verbose Print output to the console | 打印输出到控制台
# @importFrom pbapply pblapply
# @return return a list containing nn index and nn multimodal distance | 返回一个list，包含最近邻下标和多模态距离
#
#' @importFrom methods new
#' @importClassesFrom SeuratObject Neighbor
#
MultiModalNN <- function(
  object,
  query = NULL,
  modality.weight = NULL,
  modality.weight.list = NULL,
  k.nn = NULL,
  reduction.list = NULL,
  dims.list = NULL,
  knn.range = 200,
  kernel.power = 1,
  nearest.dist = NULL,
  sigma.list = NULL,
  l2.norm =  NULL,
  verbose = TRUE
){
  # (A1) 实例T，所以输出进度条
  my.lapply <- ifelse(
    test = verbose,
    yes = pblapply,
    no = lapply
  )

  #(A2) 设置参数：基本都是从 FindModalityWeights() 返回的对象中的 params slot 中获取
  k.nn <-  k.nn %||% slot(object = modality.weight, name = "params")$k.nn #20
  reduction.list <- reduction.list %||%
    slot(object = modality.weight, name = "params")$reduction.list
  #$pca
  #[1] "pca"
  #
  #$apca
  #[1] "apca"
  dims.list = dims.list %||%
    slot(object = modality.weight, name = "params")$dims.list
  #1:30, 1:20
  nearest.dist = nearest.dist %||%
    slot(object = modality.weight, name = "params")$nearest.dist #距离每个细胞最近的一个细胞的欧氏距离
  #> modality.weight_2@params$nearest.dist |> names()
  #[1] "pca"  "apca"
  sigma.list =sigma.list %||%
    slot(object = modality.weight, name = "params")$sigma.list #每个模态的 bindwidths
  l2.norm = l2.norm %||%
    slot(object = modality.weight, name = "params")$l2.norm #是否对PC做L2标准化，T
  modality.weight.value  <- modality.weight.list %||%
    slot(object = modality.weight, name = "modality.weight.list") #每个模态的权重：重点
  names(x = modality.weight.value) <- unlist(x = reduction.list) #给list设置 name

  #(A3) 如果是 Seurat 对象，则获取PCA细胞嵌入：前若干个PC
  if (inherits(x = object, what = "Seurat")) {
    reduction_embedding <- lapply(
      X = 1:length(x = reduction.list),
      FUN = function(x) {
        Embeddings(object = object, reduction = reduction.list[[x]])[, dims.list[[x]]]
      }
    )
  } else {
    reduction_embedding <- object
  }

  #(A4) 实例查询为空，则 设置为PCA降维
  if (is.null(x = query)) {
    query.reduction_embedding <- reduction_embedding
    query <- object
  } else {
    if (inherits(x = object, what = "Seurat")) {
      query.reduction_embedding <- lapply(
        X = 1:length(x = reduction.list),
        FUN = function(x) {
          Embeddings(object = query, reduction = reduction.list[[x]] )[, dims.list[[x]]]
        }
      )
    } else {
      query.reduction_embedding <- query
    }
  }

  #(A5) 需要对PC做L2标准化
  if (l2.norm) {
    query.reduction_embedding <- lapply(
      X = query.reduction_embedding,
      FUN = function(x)  L2Norm(mat = x)
    )
    reduction_embedding <- lapply(
      X = reduction_embedding,
      FUN = function(x) L2Norm(mat = x)
    )
  }

  #(A6) 第一个降维的行数：细胞总数
  query.cell.num <- nrow(x = query.reduction_embedding[[1]])
  # 降维的个数：实例是2个，pca, apca
  reduction.num <- length(x = query.reduction_embedding)

  #(A7) 消息：找多模态最近邻
  if (verbose) {
    message("Finding multimodal neighbors")
  }

  #(A8) 每个模态内，找200个最近邻
  reduction_nn <- my.lapply(
    X = 1:reduction.num, #遍历模态
    FUN = function(x) {
      nn_x <- NNHelper(
        data = reduction_embedding[[x]],
        query = query.reduction_embedding[[x]],
        k = knn.range, #默认找200个最近邻？
        method = 'annoy',
        metric = "euclidean"
      )
      return (nn_x)
    }
  )

  #(A9) 获取每个模态内的200个最近邻的下标，去掉第一列(它自己)
  # union of rna and adt nn, remove itself from neighobors
  reduction_nn <- lapply(
    X = reduction_nn,
    FUN = function(x)  Indices(object = x)[, -1]
  )

  #(A10) 对每个细胞，求两个模态的200个最近邻细胞下标的并集
  # list长度就是细胞个数，每个长度不确定：最少200个，最多400个
  nn_idx <- lapply(
    X = 1:query.cell.num , #按细胞遍历
    FUN = function(x) { #x是细胞编号
      Reduce(
        f = union,
        x = lapply( #
          X = reduction_nn, #内循环遍历每个模态，y
          FUN = function(y) y[x, ] #返回这一行的细胞下标
        )
      )
    }
  )
  #求并集:
  #> union(c(1,2,3), c(3,5))
  #[1] 1 2 3 5


  #(A11) 计算细胞到其所有多模态最近邻的cell的欧氏距离：两细胞的欧氏距离 - 最近一个细胞的距离，取ReLu
  # calculate euclidean distance of all neighbors
  nn_dist <- my.lapply(
    X = 1:reduction.num,
    FUN = function(r) {
      nndist <- NNdist(
        nn.idx = nn_idx, #两模态的最近邻下标的并集，list格式
        embeddings = reduction_embedding[[r]], #原自己PC空间L2标准化后
        query.embeddings = query.reduction_embedding[[r]], #原自己PC空间L2标准化后
        nearest.dist = nearest.dist[[r]] #到最近的一个细胞的距离
      )
      return(nndist)
    }
  )

  #(A12) 模态权重距离，就是 bandwidths
  # modality weighted distance
  # 如果模态1的长度为1：实例不是长度1
  if (length(x = sigma.list[[1]]) == 1) {
    # 重复多次：和obj等长，列数/细胞个数
    sigma.list <- lapply(X = sigma.list, FUN = function(x) rep(x = x, ncol(x = object)))
  }

  #(A13) 求 加权最近邻相似性：normPC空间的距离 * 模态权重。 重要的一步，体现了用每个细胞的距离和多模态权重 计算相似性。
  # 双模态中最近邻各200个细胞（总数200~400）： exp( - 相对距离 / 模态宽度) * 模态权重
  nn_weighted_dist <- lapply(
    X = 1:reduction.num, #对模态遍历，r=模态
    FUN = function(r) {
      lapply(
        X = 1:query.cell.num, #对细胞遍历，x=cell id
        FUN = function(x) {
          # kernel.power=1
          # 该模态内(r)，该细胞的(x): exp(- 最近邻细胞相对距离 / 模态宽度) * 模态权重
          exp(-1*(nn_dist[[r]][[x]] / sigma.list[[r]][x] ) ** kernel.power) * modality.weight.value[[r]][x]
          # 分析该 加权最近邻距离：
          # nn_dist是原始相对距离，值d越小，带入指数函数 exp( - d/width) 后整体值越大
          # 后半部分 权重越大，最后整体值越大
          # 两者相乘，得到 加权最近邻相似性，该值越大，表示细胞之间的联系越紧密：细胞距离近，该模态权重高 至少一个成立
        }
      )
    }
  )
  #每个模态一个list，每个模态list长度是细胞个数，内容length不一定


  #(A14) 把各个模态的 加权相似性，按细胞相加
  # 每行长度不一
  nn_weighted_dist <- sapply(
    X = 1:query.cell.num, #对细胞遍历，x=cell id
    FUN =  function(x) {
      Reduce(
        f = "+",
        x = lapply(
          X = 1:reduction.num, #对模态遍历，r=模态
          FUN = function(r) nn_weighted_dist[[r]][[x]] #取出每个模态内，该cell的 权重距离，按对应细胞内部相加
        )
      )
    }
  )
  #> iris[1:2, 1:4] #假设输入是2行4列数据框
  #  Sepal.Length Sepal.Width Petal.Length Petal.Width
  #1          5.1         3.5          1.4         0.2
  #2          4.9         3.0          1.4         0.2
  #
  #> Reduce(f="+", x=lapply(X=1:2, FUN=function(x){
  #   iris[x, 1:4]
  # })) #
  #  Sepal.Length Sepal.Width Petal.Length Petal.Width
  #1           10         6.5          2.8         0.4


  #(A15) 权重距离 遍历细胞，按距离倒排序，给出元素下标，最大的在最前
  # select k nearest joint neighbors
  select_order <- lapply(
    X = nn_weighted_dist, #遍历细胞
    FUN = function(dist) {
      order(dist, decreasing = TRUE)
  })

  #(A16) 遍历细胞：取出 加权最近邻距离 最大的20个细胞的 下标
  select_nn <- t(x = sapply(
    X = 1:query.cell.num,
    # 下标中，取出某一行；
    # 按降序排列：加权距离的降序
    # 取出前k.nn=20个细胞的下标
    FUN = function(x) nn_idx[[x]][select_order[[x]]][1:k.nn]
  ))
  #> select_nn[1:2]
  #[1] 2014 3135

  #(A17) 遍历细胞：取出 加权最近邻距离 最大的20个细胞的 距离
  select_dist <- t(x = sapply(
    X = 1:query.cell.num,
    FUN = function(x) nn_weighted_dist[[x]][select_order[[x]]][1:k.nn])
  )
  # > select_dist[1:2]
  #[1] 0.9101496 0.8260576
  # select_dist |> hist(n=100) #[0.4, 1] 之间，正态分布，主峰在0.7

  #(A18) 换算成更容易懂的距离：越大距离越远。 换算前的应该叫做 多模态相似性。
  # 1-x 再除以2。
  # 小于0的记作0，大于1的记作1。
  # 开方
  select_dist <- sqrt(x = MinMax(data = (1 - select_dist) / 2, min = 0, max = 1))
  # select_dist |> hist(n=100) #[0.2, 0.55] 之间，正态分布，主峰在0.4

  #(A19) 返回值 邻居对象
  weighted.nn <- new(
    Class = 'Neighbor',
    nn.idx = select_nn, #20个多模态最近邻cell 下标
    nn.dist = select_dist, #20个多模态最近邻cell 距离
    alg.info = list(),
    cell.names = Cells(x = query) #cell id
  )

  return(weighted.nn)
}








(13) NNdist()：求每个细胞到最近邻的细胞的距离
$ find . | grep -v "dustbin/" | grep -v "notes/" | xargs grep -n "NNdist" --color=auto 2>/dev/null
./seurat-4.4.0/R/clustering.R:1556:NNdist <- function(


# Calculate NN distance for the given nn.idx
# @param nn.idx The nearest neighbors position index
# @param embeddings cell embeddings
# @param metric distance metric
# @param query.embeddings query cell embeddings
# @param nearest.dist The list of distance to the nearest neighbors
#
NNdist <- function(
  nn.idx,
  embeddings,
  metric = "euclidean",
  query.embeddings = NULL,
  nearest.dist = NULL
) {
  #(A1) 是list
  if (!is.list(x = nn.idx)) {
    nn.idx <- lapply(X = 1:nrow(x = nn.idx), FUN = function(x) nn.idx[x, ])
  }
  #(A2) 传入的有
  query.embeddings <- query.embeddings %||% embeddings

  #(A3) 快速计算x和y的距离，每个细胞最近邻的下标并集 n=nn.idx?
  nn.dist <- fast_dist(
    x = query.embeddings,
    y = embeddings,
    n = nn.idx
  )

  #(A4) 非空
  if (!is.null(x = nearest.dist)) {
    nn.dist <- lapply(
      X = 1:nrow(x = query.embeddings), #遍历细胞
      FUN = function(x) {
        r_dist = nn.dist[[x]] - nearest.dist[x] #减去最近一个细胞的距离
        r_dist[r_dist < 0] <- 0 #按ReLu返回
        return(r_dist)
      }
    )
  }

  return(nn.dist)
}









(14) fast_dist() C++函数
$ find . | grep -v "dustbin/" | grep -v "notes/" | xargs grep -n "fast_dist" --color=auto 2>/dev/null
./seurat-4.4.0/src/fast_NN_dist.cpp:34:List fast_dist(NumericMatrix x, NumericMatrix y, List n) {


// [[Rcpp::export]]
List fast_dist(NumericMatrix x, NumericMatrix y, List n) {
  // extracting the number of element in the knn graph
  size_t ngraph_size = n.size(); #n的长度

  if (x.nrow() != ngraph_size) { return List(); } #如果长度不等，返回空list

  List distances_list = clone(n); #克隆n这个列表
  
  // looping over the neigbors #遍历每个最近邻的key，就是 cell id
  for (size_t i=0; i<ngraph_size; ++i) {
    // extracting the neighbors vector
    SEXP ll = n[i]; #每个细胞的 双模态的最近邻细胞的下标
    NumericVector neighbors(ll); #声明向量：细胞下标，长度在200-400之间
    
    // prepopulating the distance vector #声明距离向量：和邻居等长
    NumericVector distances(neighbors.size());
    
    // extracting the row iterator of x
    NumericMatrix::Row row1 = x.row(i); #x的第i行：该细胞的每个normPC的值

    for (size_t j=0; j<neighbors.size(); ++j) { #内循环，遍历下标
      size_t n_idx = neighbors[j] - 1; #R中1-based，-1得到0-based
      
      // extracting the row iterator of y
      NumericMatrix::Row row2 = y.row(n_idx); #y的第n_idx行：该细胞的每个normPC的值
      // extracting the distance
      double distance = euclidean_distance(row1.begin(),
                                           row1.end(),
                                           row2.begin()); #计算欧氏距离，实现见(15)
      
      if (distance == -1) { return List(); } #如果是-1，返回空list
      distances[j] = distance; #记录到 距离向量
    }
    
    // updating the distance
    distances_list[i] = distances; #更新细胞i的最近邻cell的距离向量
  }

  
  return distances_list;
}








(15) euclidean_distance() C++函数： C++计算向量的欧氏距离
#include <Rcpp.h>
using namespace Rcpp;

// code in-parts taken from https://gallery.rcpp.org/articles/parallel-distance-matrix/
# 并行的距离矩阵

// Assumption: the end position of vector2 is implied by the end position of vector1
# 假设 vector2 的结束位置是根据vector1推断的，假设它们等长。

// generic function for euclidean distance #泛型函数

template <typename InputIterator1, typename InputIterator2>
inline double euclidean_distance(InputIterator1 begin1, 
                                 InputIterator1 end1, 
                                 InputIterator2 begin2) {
  // value to return
  double rval = 0; #距离初始化0
  
  // set iterators to beginning of ranges #迭代器变量，指向起始位置
  InputIterator1 it1 = begin1;
  InputIterator2 it2 = begin2;
  
  // for each input item #遍历每个元素
  while (it1 != end1) {
    
    // take the value and increment the iterator #取值，并且迭代器递增
    double d1 = *it1++; 
    double d2 = *it2++;
    
    // update the distance #距离递增：差的平方和
    rval += pow(d1-d2, 2); 
  }
  
  return sqrt(rval); //#开方
}



